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Non-Markovian processes are widespread in natural and lruman-nrade systems, yet explicit model¬ 
ling and analysis of such systems is underdeveloped. We consider a non-Markovian dynamic network 
with random link activation and deletion (RLAD) and heavy tailed Mittag-Leffler distribution for the 
inter-event times. We derive an analytically and computationally tractable system of Kolmogorov- 
like forward equations utilising the Caputo derivative for the probability of having a given number of 
active links in the network and solve them. Simulations for the RLAD are also studied for power-law 
inter-event times and we show excellent agreement with the Mittag-Leffler model. This agreement 
holds even when the RLAD network dynamics is coupled with the susceptible-infected-susceptible 
(SIS) spreading dynamics. Thus, the analytically solvable Mittag-Leffler model provides an excel¬ 
lent approximation to the case when the network dynamics is characterised by power-law distributed 
inter-event times. We further discuss possible generalizations of our result. 


I. INTRODUCTION 

Non-Poisson temporal statistics where time intervals 
between isolated, consecutive actions are typically not ex¬ 
ponentially distributed, seem to be the norm rather than 
the exception for many systems: For example, period of 
infectiousness [1], inter-order and inter-trade durations 
in financial markets [2], socio-networks, including emails 
[3, 4], phone calls [5], or individual-to-individuals con¬ 
tacts being fluid [6, 7]. The absence of the robust tools 
and mathematical machinery of Markovian theory is the 
source of many challenges in modelling and analysis of 
non-Markovian systems. The burst in research activity 
that successfully combines networks and non-Markovian 
processes stems from the need to develop more realistic 
models and new analytical tools. Notable examples in¬ 
clude studying non-Poisson dynamics of networks [8] and 
non-Markovian epidemics on networks [9-11]. 

The non-Markovian property is particularly pervas¬ 
ive when considering the dynamics of time-evolving net¬ 
works, be it with fast or slow timescale [12, 13]. Deriving 
simple, solvable paradigm models can facilitate progress 
in developing new mathematical tools and methods for 
analysis and increases our understanding of the true im¬ 
plications of non-Markovianity for complex systems. Em¬ 
pirically, it turns out that many inter-event distributions 
have power-law tails (see [14] and references therein). 
Therefore, it is also necessary to develop methods able 
to deal with such distributions. 

It is now widely accepted that human contact patterns 
are highly dynamic and may evolve concurrently with an 
epidemic; many Markovian models for this setup exists 
[15-17]. Here, we take the next step and consider a dy¬ 
namic network with non-exponential waiting times with 
consecutive updates which are either link activation and 
deletion [17]. As a first step in the rigorous analysis of 
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networks with Non-Markovian dynamics, we consider a 
random link activation-deletion (RLAD) model that nat¬ 
urally leads to a stochastically evolving network [17, 18]. 
This model amounts to considering undirected and un¬ 
weighted networks, where an event consists of selecting a 
link at random, independently of whether present or not, 
followed by its activation, if the link is absent, or dele¬ 
tion if the link is active. Such operations are separated 
by inter-event times sampled from the Mittag-Leffler dis¬ 
tribution, that allows for analytical tractability. This ex¬ 
actly solvable model of non-Markovian network dynamics 
is an important special case of a more general theory for 
non-Markovian processes outlined in [19], and it is re¬ 
lated to recent outstanding developments in probability 
theory [20, 21], Indeed, we provide a bottom-up deriva¬ 
tion for the master equation of some fractional birth and 
death processes in a finite capacity system, introduced in 
[20]. This allows us to compute theoretically the exact 
distribution of the total number of links in the network 
at any time and its large-time limit. We demonstrate the 
power of the analytical model by comparing it with sim¬ 
ulations using more widely-used power-law distributed 
times. The rigorous analysis of this model, including ex¬ 
plicit expressions for the distribution of the number of 
links in the network for t > 0, is followed by consider¬ 
ing a Markovian SIS epidemic on our non-Markovian 
dynamic network. Finally, we briefly discuss the gener¬ 
alization of our method to general Markov chains with 
random state changes occurring according to a generic 
renewal process. 


II. AN EXACTLY SOLVABLE MODEL 
A. Basic ingredients 

Consider an arbitrary graph on N nodes as an initial 
state of the dynamics. We are interested in the number of 
(unique undirected) links in the network at a given time 
t. We denote this number by X(t), and it takes values in 
S = {0,1,..., M} where M = N(N — l)/2, the maximal 
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possible number of links. The time periods where X(t) 
remains constant are called sojourn times or inter-event 
times. We assume that sojourn times {Tj}i>i are drawn 
independently from the family of Mittag-Leffler distribu¬ 
tions with parameter (or order) /J £ (0,1) [22]. Their 
cumulative distribution function (c.d.f.) is indexed by 
this j3 and it is given by 

F^\t) =¥{T< t} = 1 -E p (-t p ). (1) 

Here Ep{z) is the Mittag-Leffler function, defined by 


Mittag-Leffler sojourn times lead to a simpler analyt¬ 
ical treatment of non-Markovianity in the presence of 
extreme power-law tails than its cognate Pareto distri¬ 
bution. However, we do explain below how the theoret¬ 
ical framework developed here can be used to approxim¬ 
ate the behaviour of non exactly-solvable systems with 
Pareto power-law distribution, as it is most commonly 
used. For this we must introduce a scaling parameter 
(time change) 7 > 0 for the waiting times: We say that 
a random variable T is Mittag-Leffler 7 (/J) distributed if 
and only if 


E 0 {z) = 


y _: 

^ r(i 

n =0 v 


■Pn)' 


( 2 ) 


Ep is entire for all /3 > 0. At /3 = 0 the series converges 
uniformly only on a disc of radius 1 , though the function 
can be extended analytically onC\{l}. Equations (1), 
(2) define a proper c.d.f. only when /3 £ (0,1]. This 
is equivalent to the claim that, for /? £ ( 0 , 1 ], Ep{—t^) 
is completely monotone. A 00 ) function f(t) is 

completely monotone if (—1 ) n d n f(t)/dt n > 0 for all non¬ 
negative integer n and all t > 0. Now, Mainardi and 
Gorenflo [23] proved that, for /3 £ (0,1), Ep{—t^) can be 
written as a mixture of exponential distributions given 
that 


where 


/*OC 

Epi-t 13 ) = / exp (~rt)Kp(r)dr, 
Jo 


Kp{r) = 


1 r 13 1 sin(/ 37 r) 

7 t r 2 P + 2 r 13 cos(/?7r) + 1 ’ 


( 3 ) 

( 4 ) 


and 


Kp(r) dr = 1. (5) 

Therefore, complete monotonicity of Ep(—t ;3 ) is an im¬ 
mediate corollary of Bernstein’s theorem [24, 25] . A dir¬ 
ect proof that Ep(—x) is completely monotone can be 
found in reference [26]. When 0 < /? < 1 these distribu¬ 
tions are heavy-tailed with infinite mean while at j3 = 1, 
T is mean 1, exponentially distributed. This family of 
distributions interpolates between a stretched exponen¬ 
tial for small t and a power-law for large t [23]. Namely, 
one has 



Ep(-fP) ~ exp(-t / 3 /T(l + /?)), t < 1 , 

n / j.0\ sin(/ 37 r) T(/3) ^ 

Ep{~t p ) --- jp~, t^oo. ( 6 ) 


Therefore, the use of these distributions is more general 
than it might seem at a first glance. A word of notational 
caution: Here /3 is the order of the polynomial decay 
of the survival function, but most commonly power-law 
distributions are identified by the order of decay of their 
densities, which in our case is 1 + /3 £ ( 1 , 2 ). 


F^\t) = F{T < t} = 1 - Epyt/^). (7) 

For 7 = 1 the c.d.f. is reduced to that of equation (1) 
and we see that T is Mittag-Leffleri(/3) if and only if 7 T 
is Mittag-Leffler 7 (/3). 

For the rigorous derivation of the evolution equations, 
we restrict for clarity to the 7=1 case and remark 
how the equations behave with the extra scaling later. 
Fix a parameter /3 £ (0,1). The network evolves in a 
semi-Markov way: Let T\, T 2 ,... be independent Mittag- 
Leffler (/3) times and define the partial sum 

n 

S n = Y J T U n> 1. (8) 

k =1 


The sequence Si, S 2 , ■ ■ ■ denotes the event times at which 
the state of the network X(t) attempts to change. A 
change in the state means an undirected link is either 
deleted or activated. For extra flexibility, the model is 
introduced with an extra delay parameter a £ [ 0 , 1 ), so 
that if a ^ 0 allows the active links to remain unchanged 
even if there is an attempt of a change. 

It is useful to define the embedded Markov chain for 
the number of links in the network, X n ,n > 1, with 
state space S. Initially Xq = i , as we start with i present 
links and the number of links in the network increases, 
remains or decreases according to the following transition 
probabilities 


qk,k -1 = Po{X j+ i = k - l\Xj = fc} 

0 , k = 0 , 

= < 1 - a, k = M, (9) 
(1 — cx)jj, otherwise, 


Qk,k = a, 


and 


( 10 ) 


Qk,k +1 — Po{Xj+i — k + 1 | Xj — k} 

1 — a, k = 0 , 

= { 0, k = M, 

1 — a — , otherwise. 


( 11 ) 
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In words, at the time of the z-th event, we pick a link uni¬ 
formly at random out of all available links. With prob¬ 
ability a nothing changes, otherwise on the event that a 
change will happen in the system, we delete or add a link 
in the following way: If the link was active (present) in 
the network, it is now deleted, otherwise it is now activ¬ 
ated. Notice that the embedded dynamics are equivalent 
to the a-delayed version of the Ehrenfest chain. 

To connect the embedded chain X n with process X ( t ), 
define the counting process 

Np(t) = maxjn € N : S n < t} (12) 

that gives the number of events up to a finite time horizon 
t. This process is also called a fractional Poisson process. 
Then we have 


X{t) = X Np{t) =X n l{S n <t<S n+1 }, (13) 

i.e. the state of the process at time t is the same as that 
of the embedded chain after the last event before time t 
occurred. 


B. Semi-Markov Master Equation 

All information about X(t) is encoded in the pairs 
{( X n , T n )} n >i which are a discrete-time Markov renewal 
process, satisfying 

P{.Y n+1 = j,T n+1 < t\{X 0 , S 0 ), ...,(X n = i, 5„)} 

= lP{^n+i = j, T n +i < t\X n = i}. (14) 

Af(-) is then a semi-Markov process subordinated to 
Np(t) [18] and satisfies the forward equations 

Pi,j ( t ) = Ft > ^ gej [ Pi,i{u)f^\t - u) du. (15) 

tes 

Incidentally, a semi-Markov process is Markovian if and 
only if the distribution of {T n } n >i is exponential [27]. 
Above we introduced Pij{t) = P{Af(i) = j|X(0) = *}, 
the tail (complementary cumulative distribution func¬ 
tion) (t) = 1 — Fj?\t) and the Mittag-Leffler 

density or order /3. These equations are proved by con¬ 
ditioning on the time of the last event before time t. By 
taking Laplace transforms in (15), and using the known 
Laplace transform of the Mittag-Leffler survival function 
and probability density function 

C (EtV);s) = and C 

(16) 

followed by some straightforward algebra, the evolution 
equations for Pij(t) become (see Appendix A), 

= (17) 

+ l 1 - ») ( M A ) +l '" .' + ^PUti(*)) • 


Similarly, the equations of the boundary terms are 

d Pi${t) = _ a ) ^- Pii0 (t) + jjPi, 1 (t)') (18) 


d p p iy M{t) 

dtP 


= (1 - a) ( -p ilM (t) + jjPi,m- i(*) 


(19) 


Symbol /dt& in (17), (18), (19), denotes the /3 frac¬ 
tional Caputo derivative [22] of a function fit) given by 


d 0 f(t) 

dtfi 


1 


T(l-/3) 




When (0 = 1, equations (17), (18), (19) reduce (as ex¬ 
pected) to the standard Kolmogorov equations for the 
Markovian RLAD [17]. These equations also explain 
analytically why a is called the delay parameter. When 
considering the scaled Mittag-Leffier 7 (/3) times, equation 
(17) becomes 


d^Pi^j (t) 

dtP 

+ 7" ,9 (1 


-7 ^(l -a)p^ht) 


t,3 


- a) 


M-j + 1 

w 


Pu-i(t) 


( 20 ) 


.7+1 

M 



and similarly for the boundary equations. Specifically we 
see, as in the Markovian case, that a scaled sojourn time 
distribution results in a (fractional) scalar multiple of the 
forward equations. 


C. Exact solution 

Equation (15) gives an analytical way to obtain the 
fractional equation for the evolution of the transition 
probabilities, but it is not very useful for computational 
purposes. Instead, it is fruitful to find the solution of the 
system of equations (17), (18), (19) by a simple condi¬ 
tioning argument on the values of Np(t) (see Appendix 

B) 

OO 

PiA*)=rfPm+Q$nmt) = *}. m) 

n—1 

where q^ are the n-step transitions of the embedded 
discrete Markov chain, namely the entries of the n-th 
power of the transition matrix Q defined by equations 
(9), (10) and (11). The distribution of the fractional 
Poisson process has a simple expression generalising the 
Poisson distribution [28], namely 

fP n t \ 

V{NAt)=n} = —E$ l \-tP), (22) 

where E^{—tP) denotes the n-th derivative of Ep{z) 
computed for z = — 1@ . Equation (21) can also be verified 
to satisfy (17) using Laplace transforms, (see Appendix 
B). 
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III. RESULTS AND APPLICATIONS 

All simulations are event driven, both for dynamic net¬ 
works and when this is coupled with epidemic dynamics. 
Waiting times for all the possible events are generated 
from appropriate distributions. Hence, the next change 
or an update is always determined by the smallest waiting 
time and the event corresponding to it is executed. This 
is then followed by the necessary update of the waiting 
times of the events affected by the most recent change. 
In reference [29] , readers can find an alternative efficient 
simulation method which effectively extends the ideas of 
the Gillespie algorithm from the Markovian to the non- 
Mar kovian case. 


A. Explicit calculation of Pi : j(t) 

The probabilities involving the counting process Np ( t ) 
have an explicit integral representation [30] and for nu¬ 
merical purposes they can be approximated well either 
with Monte Carlo simulations or with a numerical in¬ 
tegration scheme. Once the transition probabilities of 
the embedded Markov chain are known, every term is 
known in (21) and it can be used to exactly calculate the 
non-equilibrium probabilities Pij(t) (see Appendix B). 
The excellent agreement between theory and simulation 
is shown in Fig. 1. 



Link Count 


Figure 1. (Color online) Comparison between Monte Carlo 
simulations and theory. The discrete markers are the es¬ 
timated probabilities pi9o,j(250), averaged over 10000 Monte 
Carlo simulations starting from a fully connected network 
with N = 20 nodes and for /? = 1, 0.7, 0.5, as we move from 
left to right. The Monte Carlo simulations were performed us¬ 
ing an event-driven algorithm taking non-Markovianity into 
account. The solid curves are the theoretical predictions as 
dictated by equation (21). 

An immediate application is to use equation (21) 
and compute theoretically and numerically the expected 
number of active links in the network at a given time. 
Starting from any initial number of active links *o, use 


( 21 ) to compute 

E l °(X(t)) =E'°(elQ N ^v 0M ) 

= ef o r °{Q n ^)v 0 ,m- (23) 

In the equation above e ; 0 is the standard basis factor with 
1 at the io-th coordinate and v 0 ,m = (0,1,2,..., M) T . 
Note that in the particular case where Q diagonalizes, the 
analytical expression for the expectation (23) is merely 
a linear combination of different values of the probabil¬ 
ity generating function of Np(t), G^\s;t), given by (see 
[31]) G^\z-t) = E (z N ^) = Ep((z - 1)V). Note that 
when Q diagonalises, there is no need for simulating a 
large number of realisations to estimate the expectation; 
a fast numerical integration scheme is sufficient. 

B. Approximation of Pareto-distributed 
inter-event times 

We now compare the behaviour between two RLAD 
networks; one with Mittag-Leffler times and one where 
we alter the waiting time distribution to a generalised 
Pareto(^) with density 

Mt) = t > 0. (24) 

Exponent S = 1 + /3 in order for the tails of Mittag-Leffler 
and the Pareto to have the same behaviour at infinity. 
In fact we compare the two networks over three layers 
of increasing complexity. First, in Fig. 2(b,c) we plot 
the probability mass function for the number of singly- 
counted links at a pre-specified time horizon T = 2000, 
averaged over 5000 simulations. As a point of reference, 
output from the Markovian RLAD is shown in Fig. 2(a). 
The theoretical curve at equilibrium is the large time 
limit of equation ( 21 ) and it is the mass function of a bi¬ 
nomial distibution with M trials and success probability 
1/2. Second, in Fig. 2(d) we plot E(A(f)) as a function 
of time up to time T = 2000. The two curves could also 
be computed based on equation (23). 

The excellent agreement is achieved by finding a suit¬ 
able scaling 7 so that the c.c.d.f. (survival functions) of 
the two distributions match well, at least up to the pre¬ 
specified time horizon that we want to study. Further 
details can be found in Appendix C. The matching is 
good for (3 < 0.9 while for larger /3, this idea can be 
used to study the stochastic dominance between the two 
coupled networks and offer rigorous bounds. 

C. Markovian SIS on non-Markovian RLAD 

Finally, we compare the two network dynamics indir¬ 
ectly, when we allow a Markovian epidemic to run while 
the networks evolve. As discussed above, human activity 
tends to be bursty and non-Markovian [32]. During an 
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epidemic, individuals become wary of the risk posed by 
it and one way to avoid infection is by limiting or redu¬ 
cing their number of contacts. This justifies the deletion 
of links as time evolves. On the other hand, close con¬ 
tacts cannot realistically be removed and some level of 
communication and social cohesion must be maintained. 
Such behaviour in activation-deletion is not necessarily 
Markovian in nature, thus alternative non-Markovian dy¬ 
namic network models are necessary. 

Nodes in the network represent individuals from a pop¬ 
ulation and links describe the contact patterns amongst 
these. Each individual can be either infected (/) os 
susceptible ( S ). An infected individual remains infec¬ 
ted for exponentially distributed periods of time Th he. 
Th ~ Exp(l/r#), where th is the average time in which 
infectious individuals are healed. Similarly, infection oc¬ 
curs at the points of a Poisson process with time to infec¬ 
tion Tj exponentially distributed, i.e. Tj ~ Exp(l/r/), 
where tj is the average time in which an infection spreads 
across a link connecting a susceptible and an infected 
node. In this framework, both network and epidemic dy¬ 
namics can be considered in the context of event-driven 
simulations, where the timing of the next state change 
is always determined by the smallest waiting time and 
the precise event corresponding to it. The epidemic does 
not interfere with the network dynamics, however its 
propagation is intertwined with the background dynamic 
network topology. Initially, before the infection starts 
spreading, we assume that all links are present, in order 
to avoid early stochastic extinction. 

The simulations in Fig. 2 show the prevalence (pro¬ 
portion of infected individuals (Fig. 2(e)) on a Mittag- 
Leffier 7 (/3) RLAD network (solid lines) and a direct com¬ 
parison (square markers) with the Pareto(d). Again, we 
use the same sets of /3,y ,5 as before and we emphasise 
the excellent agreement between the two. 

Incidentally, as expected, the non-Markovian network 
dynamics create a striking effect by slowing down the 
network dynamics and thus effectively blocking the at¬ 
tainment of statistical equilibrium in a realistic time ho¬ 
rizon (Fig. 1, 2(a,b,c)). This leads to a heightened level 
of infectiousness in the population (Fig. 2(e)) and high¬ 
lights the importance of quick reactions. Naturally the 
statistical equilibrium will be reached after a much longer 
time period, but the delayed curves can now be theoret¬ 
ically computed or approximated. One way to explain 
this delayed convergence to equilibrium is to look at the 
mixing time of the embedded chain in the total variation 
distance. To be more specific, the number of active links 
is a continuous time, irreducible birth-death chain with 
a unique binomial invariant distribution i r, independent 
of the delay parameter a, given by 

7r fc = lim P{I(t) = k\X(0) =i}= ( M ) 2 -M , k £ 5. 
t-v oo \ k J 

(25) 

This can also be deduced from the fact that in the aperi¬ 
odic a-delayed case, at equilibrium, individual graphs are 
uniformly distributed. That is because the chain on the 


set of distinct graphs has a doubly stochastic transition 
matrix. With this in mind, the degree distribution of 
a single node in network chosen uniformly at random 
can be immediately computed as follows. Let iq be a 
selected node in G and define G Vl the subgraph of G 
where v\ and all its incident links are deleted. G Vl is 
now a graph on the set {v 2 ,... ,wjv} that has at most 
K = ( Ar 2 ~ 1 ) = M — N + 1 links. Let deg(rq) denote the 
degree of tq. Then, 


Pjdeg^i) = £} = ^ 2 M 

graphs G: deg ( vi)=£ 

= 2 _M card{graphs G : deg (iq) = £} 
r N - P 


= 2 


-M 


card{graphs G Vl } 


_ 2~k~n+i 


N — 1 

£ 


N - 1 

£ 

2 ~ n + i. 


K 

E 

2—0 


K 


(26) 


The third equality is a counting argument. The number 
of graphs such that iq has exactly £ incident links, is con¬ 
structed by first selecting where those links go, and then 
constructing the subgraph G Vl . This can be understood 
heuristically as follows. Select a link and focus on one of 
the two nodes. If this node has h active links, this number 
will either go to h+1 with probability (N — 1 — h)/(N— 1) 
or to h — 1 with probability hj (N — 1). This leads to an 
invariant binomial degree distribution (26), with an av¬ 
erage degree of (N — l)/2, which amounts to N(N — l)/4 
active edges in the network in line with the average of 
the link distribution from (25). The chain mixing time 
tmix(e) is the minimal time so that the total variation 
distance between the measures n and p(t) is smaller than 
some tolerance e, i.e. 


Ib(tmix(e)) - tt||tv = sup |pfc(imix(e)) - 7Tfc| < e. (27) 
fees 

For the Markovian RLAD with a > 0, use Theorem 1.1 
and Example 4.3 in [33] to see 

tmix(e) < Ce~ 2 M 2 \ogM. (28) 


Thus, the Markov chain approximates relatively well its 
equilibrium by time Ce~ 2 M 2 log M. In particular this 
implies that, on average, the Markovian RLAD continu¬ 
ous chain needs 0(M 2 log M) time, and therefore by the 
law of large numbers, it needs this order many events 
until it is well mixed. In fact this bound is also true for 
the embedded discrete chain. This n = M 2 log M should 
be considered as a necessary lower bound of steps so that 
the sample average of the probabilities of the embedded 
chain approximates n. Therefore in order to have an ac¬ 
ceptable level of accuracy for the embedded chain when 
the RLAD has Mittag-Leffler waiting times, using (22), 
we need a higher polynomial order 0{M 2 ^ (log M) 1 ^) 
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Figure 2. (Color online) Top: Distribution of links at t = 2000 and comparison between Mittag-Leffler and Pareto. Panels 
(a),(b),(c) show the distribution of singly counted link numbers at t = 2000, with theoretical equilibrium prediction (continuous 
line) and simulation results (discrete markers). Theoretical equilibrium is the same for all values of /3 € (0,1] and is given by 
the binomial distribution (25) here drawn as a continuous line for ease of representation. □ markers in the figures (a), (b), (c) 
correspond to (/3,7) = (1,1), (0.7,3.14), (0.5,4), respectively and <0 markers for the corresponding Pareto(<5). Bottom: On the 
second row, we show (d) the expected number of singly counted links in the network (the dashed line is the theoretical prediction 
of N(N — l)/4 = 95 and (e) the expected prevalence. The continuous (noisy) line is for (/3,7) = (1,1), (0.7, 3.14), (0.5, 4), 
respectively from bottom to top and □ markers are the corresponding values for the matched Pareto network. The networks 
have N = 20 nodes, and simulations start with a completely connected network. The SIS epidemic is simulated as a Markovian 
process with transmission and recovery rate tj = 0.25 and th = 1, respectively. The spreading process initialises with 5 
infectious nodes. Since we start with a fully connected network and a slow network dynamics, the prevalence rapidly increases 
from 25% to almost 80%. This also reflects the relation between the time scales of the network and epidemic dynamics, with 
the epidemics being much faster in this example. The simulation is event-driven and is implemented by keeping track of all 
events and their waiting times. These are averaged over 5000 simulations and use a = 0 (periodic case), so all events create a 
change in the network and equilibrium is reached earlier. 


time in order to guarantee on average the same number 
of events, and thus to guarantee a near equilibrium be¬ 
haviour for the embedded chain. Study of the slow-down 
phenomenon for non-Markovian dynamic networks, using 
the total variation distance, can be found for example in 

[34]. 


IV. DISCUSSION AND CONCLUSIONS 

A. Generalization 

Equation (21) can be generalized to any counting 
process and any discrete Markov chain and, as a con¬ 
sequence, to any embedded Markovian graph dynamics 
[18]. To be more specific, let qij denote the one-step 
transition probability from state i to state j for a dis¬ 
crete Markov chain X n and let N(t) be a generic counting 


renewal process. Then, for the process 

X(t)=X NW =X n l{S n <t<S n+1 }, (29) 

the probabilities Pi,j(t) = P{A(f) = j|AT(0) = i} are 
given by 

OO 

Pi,j(t) = + Y Qi^{ N (t) = n}, (30) 

n =1 

where the symbols have the same meaning as in (21) 
and is a sequence of i.i.d. positive random vari¬ 

ables with the usual meaning of inter-event times with 
arbitrary distribution, not necessarily with fat tails and 
infinite mean. The reader is invited to follow the first 
proof of Appendix B by replacing Np (f) with a generic 
counting renewal process. This will convince the reader 
of the wide generality of this result. A heuristic argu¬ 
ment to justify (21) and (30) runs as follows. In the 
time interval (0, i), n > 0 events may have occurred. In 
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the case n = 0, at time t the process is still in state i 
and P {N(t) = 0} = P {T > t} = Frit). If n > 1, the 
probability of being in state j after n events is given by 
qfj ■ Given the independence between the renewal pro¬ 
cess and the Markov chain, the probability of being in 
state j at time t and n transitions occurring in the time 
interval (0,f) is q^F{N{t) = n}. Now, all these events 
are exhaustive and mutually exclusive. Then, total prob¬ 
ability and infinite additivity imply that Pij(t ) is given 
by (30). These considerations suggest further general¬ 
izations taking into account possible dependence within 
the couple {(X n ,T „)} n >i as well as serial dependence or 
state dependence of inter-event times. 


B. Example: A simple probabilistic model for 
relaxation in dielectrics 

In order to illustrate the generalization discussed above 
with an example, we consider relaxation phenomena 
[35]. Probabilistic modelling of relaxation assumes that a 
physical system (e.g. a molecule) can exist in two states 
A and B. We further assume that state A is transient and 
state B absorbing, so that the deterministic embedded 
chain has the following transition probabilities qA.A = 0, 
qA,B = 1, Qb,a = 0, and qB,B = 1- This means that if 
the system is prepared in state A, it will jump to state B 
at the first step and it will stay there forever. Now sup¬ 
pose that the inter-event time T is random and follows 
an exponential distribution with rate A = 1 for the sake 
of simplicity. Based on equation (30), we immediately 
have pa,A it) = Fx(t) = exp(—t). Therefore, the prob¬ 
ability of finding the system in the initial state decays 
exponentially towards zero. This relaxation function is 
the solution of 

4 -PA,Ait) = -PA,A(t ), Pa,a{ o) = 1. (31) 

at 

The response function is defined as £p(£) = — dpA,A(t)/dt 
and its Laplace transform is 1/(1 + s). For s = —iui this 
is the Debye model [35]. If inter-event times follow the 
Mittag-Leffler distribution, we get PA,A{t ) = i'Y(t) = 
Ep(—tP). This is the solution of [28] 

dP 

■J t pPA,A{t ) = - PA,A(t ), PA,A( 0) = I- (32) 

In this case, the Laplace transform of the response func- 
tion£cc(0 = —dpA,A(t)/dt is l/(l+s /3 ) and for s = — iu), 
we get the Cole-Cole model [35-37]. 

C. Final considerations 

In conclusion, we provide an exactly solvable non- 
Markovian dynamic network model. The R.LAD is par¬ 
ticularly attractive as it has analytical and numerical 
tractability coming from fractional calculus. We are able 


to explicitly use the master equation formalism and ana¬ 
lytically derive the distribution of the number of links in 
the network for arbitrary times X(t), consequently com¬ 
puting E(X(f)). We highlight an important connection 
and possible avenue to approximate non-Markovian prob¬ 
lems using fractional calculus, by coupling a Pareto net¬ 
work and show the agreement with the tractable model. 
Moreover, we discuss how our result can be extended to 
a generic counting renewal process. 
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APPENDIX 

In this Appendix, we cover the rigorous proofs of the 
equations shown in the main text and further clarify some 
notions. Some details about the procedure used to couple 
the Pareto distribution with the Mittag-Leffler are also 
highlighted. 

A. Derivation of fractional equations. 

We want to show that equations (17), (18) and (19) in the 
main text are obtained from (15). The analysis proceeds 
by way of Laplace transforms. They are defined as 

f-OO 

£( 9 (t);s)= dt g(t) e~ st (33) 

J o 

for a suitable function g(t). In the case of the Mittag- 
Leffler distribution defined in the main text, we have 

C (f^(£);s) = and C (/^ } (£);s) = 

(34) 

For the computation that follows we use the symbol 
g(s) to denote the Laplace transform C(g; s) of any func¬ 
tion g. Taking the Laplace transform of (15) and using 
equations (4), (5), (6) in the main text for our particular 
example, we have for 1 < j < M — 1 

~(/3) 

Pij(s) = F T (s)5ij + frf’(s)apij(s) (35) 

+ / t /3) ( s )( 1 - a ) 

(M — jTl~ / \ J + 1 . /A 

x { -+ -jj-pi, j+ i{s)j . 

The boundary cases j = 0,j = M have Laplace trans¬ 
forms 

—(.P) ~(d) (f — ^ \ 

Pi,o(s) = F t ( s)S i0 + f T (s) I — — Pi,i{ s ) +ap l)0 (s)J , 

( 36 ) 







~(/3) 

Pi,m{s) = Frp (s)SiM 

1 — a 


+ # ) 


(s) 


M 


Pi, M -i(s) + ap itM (s) (37) 


respectively. We finish the computation starting from 
(35), in the case where 1 < j < n — 1. The remaining 
cases follow similarly. Multiply both sides of (35) by s 
and then subtract Pij(O) = Sij from both sides. Then, 
using (34), equation (35) becomes 


c (^ M ; s )= s ^\ s ) Si: 


§ 

Sij - Pi,j( 0) + 1 p Oipi'jis) 
s(l - a) (M - j + 1 

+ t—-' 

„ o _ , . 

TT^ apij(s) 

+ s(l - a) ( M - j + 1 


1 -a) (M -j + \_ j + 1- , f\ 

T?~( M p,J - lis) + 

8 

l _|_^S 1 Sij _ Sij + 

l + s/3 l Af )+ M Pi ’ i+l(s) J’ 

thus, after some algebraic manipulations we have 


4%^*) _ 


-Si. 


-apij(s) 


(38) 


+ (!-“) 


M -j + 1 

w 


-Pi,j- l( s ) + ) 


Focus for the moment on the factor s 1 (1 + s^). Its 
inverse Laplace transform is 


standard law of total probability, where the space is par¬ 
titioned according to the number of jumps of the counting 
process Np(t): 

Piij (t)=¥{X(t)=j\X(0)=i} 

oo 

= Y l P{X(t)=j,N p (t) = k\X(0) = i} 

k=0 

oo 

= Y,n x (t)=j\ N p(t) = k,X(O)=i}P{N 0 (t) = k} 

k=0 

= ¥{X(t) = j\Np(t) = 0,X(0) = i}¥{N p (t) = 0} 

OO 

+ £>{*(i) = j\N p (t) = k,X( 0) = i}¥{Np(t) = k} 

k—1 

= P{X(t) = j\T x > t, X(0) = i}¥{Ti > t} 

OO 

+ ^P{X(f) = j\Np(t) = k,X( 0) = i}¥{Np(t) = k} 

k=1 

= SijF^\t) 

OO 

+ Y ¥{X k = j\l{S k < t < S k+1 }, X 0 = i} 

k=1 

x¥{N 0 {t)=k} 

OO 

= SijF^\t) + Y = 3\Xo = *} p {^/3 (t) = k}, 

k= 1 

where it finally leads to 


C 1 (s : (1 + s 13 ); t) 


t~P 

r(i — P) 


+ i = $f,{t) + i. 


(39) 


PiAt) = SijF { £\t) + Y qWF{N p (t) = k}. (41) 

k=1 


Kernel &p(t) is what is used in fractional calculus to 
define the /3 fractional Caputo derivative (see reference 
[19] in the main text) of a function /(t), given by 


dPf(t) 

dt? 





dm 

dt' 


dt’. 


Thus, use (39) to write the left hand side of (38) as 
a product of two Laplace transforms. Then take the 
Laplace inverse of (38) to conclude 


%W=-(l-o) BjW (40 ) 

+ (! - o) ( M jJ + 1 W.j-i(Q + ' 

Similarly, the equations of the boundary terms are de¬ 
rived (18), (19). 


B. Solution to the fractional equations. 


Equation (41) is equation (21) in the main text and 
as we say, gives the theoretical solution to the fractional 
equations, because of an explicit integral representation 
of ¥{Np(t) =n}. It is given by 

P{Np(t)=n} = — Ef\-lP) 

r°° / u\ u n 1 

Function Fs p ( t ; u) is the c.d.f. of a stable random variable 
Sp(u, 7 , S) with index /3, skewness parameter v = 1, scale 
7 = (u cos( 7 r/ 3 / 2)) 1 / /3 and location 5 = 0. This integral 
representation was used to numerically compute the solid 
curve in Figure 1 [30]. 

We now verify via Laplace transforms that this solu¬ 
tion (41) indeed verifies the fractional equations. For 
simplicity we set the delay parameter a = 0 and we only 
show it for equation (17). We need 

c (^l ; s ) = ~ ^ _ 1 5 (o+), 


The solution to equations (12), (13), (14) can be seen , . _ , . s (f(P), f t ( s ) 

to be equation (15) in two different ways. One is the ^ ^ /3l 1 — n r, s ) — t b s l yJr \ S )J — ^ _| 
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The Laplace transform of (17) 
s^Pijis) - s /3 ~ 1 S lJ 

= -Pij{s) + -- Pi,j- l(s) + -^-Pi,j+i(s), 

or after an algebraic manipulation 

(l+S^)Pi,j(s) — ^ ^ij l,jPi,j —1 (^)H”Q , .7+l,jPi,j+l (^)• 

(42) 

To verify (42), directly take the Laplace transform in (41) 
to write 

Pi,j(s) = SijFj. (s) + F t (3) J2 QS (/t } (s)) (43) 

fc =1 

and substitute to the right hand side in (42) that now 
reads 

+ Qj-l,jPi,j-l(s) + Qj+l,jPi,j+l{s) 

= s^~ 1 5ij 

+ <H-1,3 F T ( 8 )( 5 i,3-l+Y^Qi k j-l(fT ) ( a )) ^ 


fc=l 


Qj+h3 F T ( s ) ( kj +1 + J2 q i k j + 1 (/t 3) ( S )) 


fc=l 




— + 9i+l,j^ij' + l)^T ( s ) 

+ ^T ) (s)£(0j-i,j9i5 ) - 1 +?j+i 1 j9i5 ) +i) (/t^OO) 

~(/ 3 ) 

=1^1, -v , . \ fc + 1 


fc=l 


— + qijFrp (s) 


(i+^rwE-/. 1 ? 11 (/r<»)) 


/c—1 


= s ' 3 + q itj s p 1 fj 3 \ s ) 

OO 


+ 


\ fe +1 


(l + ,«)^^(»)E4 < ‘ +I, (/t><<)) 

/c=l 

00 , 

*« ++£«S’(/?”(«)) 


« —(-9) 

= (1 + s i3 )F t ( s ) 


fc=i 


= (l + s^)pij(s), 

which is the left hand side of (42). 


C. Stochastic coupling with the Pareto distribution. 


Now we show how the complementary cumulative dis¬ 
tribution functions of the Pareto(<5) distribution and the 
Mittag-Leffler 7 (/3) of the same exponent can be shown to 
match just by manipulating the scaling factor 7 . 

On a double logarithmic scale the c.c.d.fs have a linear 
behavior at infinity with slope 1 — 6 = —/3. Our simula¬ 
tions have a time horizon of T = 2000 so the scaling 7 is 
chosen so that the c.c.d.f’s agree well for values around 
and before the time horizon. The initial value of 7 to be 
tested for matching is the solution to the equation 

sin(/ 37 r) T(/3) 1 / n 

7 r (t/ 7 )^ f 5-1 ^ ^ \si n (/37 r )r(/3) 

that implies the agreement of the asymptotic behavior of 
the survival functions. This first 7 choice will need to be 
adjusted, depending on our choice of time horizon, but 
the match can be achieved relatively well for moderate /3 
values (/3 < 0.9) (see Figure 3). 




Figure 3. (Color online) The figure shows the c.c.d.f. of a 
Pareto distribution with tail exponent <5 and the matching 
with the corresponding Mittag-Leffler 7 (/3). The left panel 
presents the case j3 = 0.7, 7 = 4 and 8 = 1.7, whereas the 
right panel presents the case /3 = 0.5, 7 = 3.14 and 8 = 1.5. 
The Pareto distribution is drawn using diamonds whereas the 
Mittag-Leffler is drawn using a solid fine. 
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